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I. INTRODUCTION 

The properties of the quantum Heisenberg model play 
a fundamental role in the physics of many-body effects 
for models defined by quantum Hamiltonians on a lat- 
tice, in several spatial dimensions [1, 2]. One of the first 
non-perturbative methods devised to study the quantum 
Heisenberg model is known as spin- wave theory (SWT). 
This is a type of mean-field theory method that is spe- 
cially suited to study the quantum fluctuations of inter- 
acting spins. The basic assumption is the existence of a 
ground state that spontaneously breaks the global sym- 
metry of the Heisenberg Hamiltonian. In this case, it 
corresponds to rotational symmetry SO (3) about an ar- 
bitrary axis. In SWT this symmetry is broken by fixing 
a preferred axis called magnetization axis of the ground 
state, and excitations appear in the form of fluctuations 
from the fixed direction. These are the Goldstone bosons 
of this spontaneously breaking mechanism and repre- 
sent the magnon modes propagating as spin waves in 
the quantum system. However, spatial dimensionality 
is crucial in order to have a well-defined semiclassical ex- 
pansion in the parameter 1/5, where S is the total spin 
at each site of the system lattice. Namely, in a quantum 
antiferromagnet the spatial dimension of the lattice has 
to be large enough in order to sustain the assumption 
of a given order in the ground state. Otherwise, strong 
quantum fluctuations in one-dimensional lattices break 
the long range order and makes the SWT invalid. How- 
ever, many interesting systems are materials in 3D and 
SWT provides very good approximations to their observ- 
able quantities. 

SWT has been extensively developed in many aspects. 
It has become by now an standard and reference tool in 
order to have a good approximate description of quantum 
antiferromagnetic systems, whenever the validity of its 
application is justified. 

To the best of our knowledge, there is an important as- 
pect of SWT that remains vaguely explored. Namely, the 
modification of SWT in order to adapt it to describe the 
natural interaction of a quantum antiferromagnet with an 
external or surrounding thermal bath that is interacting 
with it. A typical example is provided by the phonons of 
the lattice where the quantum spins are located. This is 



a basic and fundamental problem since it entails the de- 
scription of both dynamical effects, i.e. time-dependent, 
as well as finite temperature effects outside the state of 
thermal equilibrium. 

Embedding thermal fluctuations in the dynamics of a 
system may be approached from several points of view. 
For instance, in the classical domain it is common to con- 
sider the effect of a noisy magnetic thermal field acting 
on the Heisenberg Hamiltonian [3]. However, that situa- 
tion is different from what we focus in this work, where 
the noise is described from a microscopic model based on 
thermal excitation of the surrounding environment. The 
branch of the quantum theory which deals with this kind 
of problems is the theory of open quantum systems [4-7] 
that plays a fundamental role in quantum information 
theory [8, 9]. From this point of view, the quantum mag- 
net is considered as an open system, which exchanges 
energy with its environment. 

The best method to describe an open system strongly 
depends on the explicit nature of each situation. In this 
work we have applied the Davies formalism, which is a 
suitable description of an open system weakly interact- 
ing with a large environment. One its main features is 
that it allows us to derive an evolution equation for any 
spin observable of the quantum antiferromagnet coupled 
to a generic thermal bath at a certain temperature T. 
Namely, it provides us with an equation for the evolution 
of density matrix p(t). Furthermore, as a consequence 
of how this fundamental equation is obtained, a series 
of interesting results for the enlarged SWT have being 
obtained: i/ the quantum antiferromagnet thermalizes 
towards the Gibbs state for long enough times; ii/ the 
decay rate of this thermalization process can be obtained 
in closed analytical form as a function of the lattice mo- 
mentum; iii/ the thermal bath cannot be arbitrary in 
order to ensure the convergence of any observable to its 
thermal value, but it has to belong to the class of super- 
ohmic baths with specific parameters depending on the 
quantum antiferromagnet; iv/ the staggered magnetiza- 
tion can be computed analytically and we can obtain its 
behaviour with time and temperature thereby unveiling 
the fate of the antiferromagnetic order parameter and v/ 
the thermal evolution of the magnon form factor can be 
also computed explicitly. These quantities are of physical 
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importance and observable in inelastic neutron scattering 
experiments [10-12], Raman experiments [13] etc. 

This paper is organized as follows: in Sect. II we review 
the linear spin-wave theory and establish our notation. 
Sect. Ill describes the microscopic coupling of the SWT 
Hamiltonian with a Hamitonian bath of bosonic opera- 
tors. In Sect. IV we derive the complete master equation 
for describing the evolution and thermal effects of the 
whole interacting system-bath. In Sect.V we compute 
relevant observables under the above conditions, such 
as the staggered magnetization, two-spin correlators and 
form factors. Sect. VI is devoted to conclusions. We refer 
to appendix A for expressions of the time-evolution of 
the first and second moments, and to appendix B for the 
detailed calculation involving the two-spin spatial corre- 
lation functions. 



II. SPIN WAVE THEORY FOR QUANTUM 
ANTIFERROMAGNETS 

First of all, let us briefly remind the spin wave theory 
for quantum antiferromagnets and to establish our nota- 
tion. The system consists of a lattice with a spin S on 
every vertex. The Hamiltonian contains only two-body 
terms between first neighbors according to a Heisenberg 
interaction 



h s = j S r 



(i) 



with J > for antiferromagnetism. The phenomenology 
displayed by this Hamiltonian strongly depends on the 
morphology of the lattice. Particularly if the lattice is 
bipartite (i.e. we can define two sublattices A and B in 
such a way that the first neighbours of a A belong to B 
and viceversa, see figure 1) the ground state is close to a 
staggered spin configuration known as Neel state. How- 
ever if the lattice is not bipartite (e.g. triangular lattice) 
the system becomes frustrated, and no simple configura- 
tion is found to be a ground state for the diagonal part of 
the Hamiltonian (1). For our purposes we shall consider 
a square lattice. 

The diagonalization of the Hamiltonian (1) is not an 
easy task, and no exact solutions are known for spatial 
dimensions d > 2, or for spins S > 1 in d = 1. Thus, 
approximation methods become very useful. Probably 
the most fundamental of them is based on the Holstein- 
Primakoff approximation [14, 15] and leads to the so- 
called spin wave theory [16], which is also applicable to 
ferromagnets [17]. This method rewrites the spin oper- 
ators in terms of bosonic annihilation and creation op- 
erators, a and a"*", [a, a"*"] = 1. Concretely, for r in the 
sublattice A 



2Sfs(a r a r )a r , 
2Sa r f s {a r a r ), 




FIG. 1: Arrange of a quantum antiferromagnet in the Neel 
state on a square lattice. The color of the spins denotes the 
two different sublattices, A (blue arrows) and B (orange ar- 
rows) . The true ground state is close to this staggered config- 
uration, however there is a slight disarrange in the orientation 
of the spins due to quantum fluctuations. 



and for r in the sublattice B 



S r 



2Sbtf s (btb r ), 
2Sfs{b%)b r , 
S, 



with 



fs(*)=(l-f s 



1/2 



(3) 



(4) 



By writing the Hamiltonian (1) at first order in fs(x) we 
obtain the so-called linear spin-wave theory: 



Hlsw = J [-NdS 2 + 2dS E r (4 a r + *>% 



(5) 



This approximation is valid to describe states where 
(fs( a t a r)} = (fs(b r b r )} ~ 1, and thus they also verify 



(a r a r ), (blb r ) <C 25. 



(6) 



This is the self-consistent condition characteristic of this 
mean-field theory method. 

The Hamiltonian Hlsw is quadratic in boson opera- 
tors, so in order to diagonalize it we take Fourier trans- 
form 



a r = 



b r — 



JLy e -ii" 



(7) 
(8) 



with A^4 = Nb = N/2 for a square lattice and the lattice 
wave vector takes on the following discretized values, 



k = 



27rra Aum 



N a ,b 



N 



(9) 



Then the first term of Hlsw is easy to compute given 



(2) the orthonormalization rule 



CL^CLip 



,reA,B* - 5 k,o- For 
the second one, we parameterize r r neighbour to r as 
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r' = r+r M , where f M is the unit vector in the fi direction, 
which in d = 3 and starting from the first site, can be 
(1, 0, 0), (0, 1, 0) or (0, 0, 1). Thus we obtain 



H 



LSW 



J 



where 



-NdS 2 + 2S Z k d{a{a k + b{b k ) 



(10) 



Next step is to perform a Bogoliubov transformation 
to new boson operators a k and f3 k , 

a k = cosh(0 fc )a fc - sinh(0 fc )/3£, (11) 
b k = -smh(0 k )ai + cosh(0 k )f3 k . (12) 

The function k is chosen so that the coefficient of a k /3 k 
and d k f^ k is zero: 



tanh(2# fc ) 



(13) 



With this choice the Hamiltonian of the system is diago- 
nalized 

H LSW = E° + J2 w(fc)(4"fc + PkPh)- (14) 



Here the energy dispersion relation is 

w(k) = 2js^d? - e k , 



(15) 



and Eq is a constant 



-JNS 



dS + 



N 



d 2 -a 



In summary, we have transformed the intricate Hamil- 
tonian (1) with interaction terms into another approxi- 
mate Hamiltonian which is just a collection of uncoupled 
harmonic oscillators, and so it is easy to write the whole 
spectrum analytically. The excitation of these harmonic 
oscillators are called "magnons" , because they represent 
the minimal collective magnetic excitation of the spin 
lattice. 



III. INTERACTION WITH A BOSONIC 
ENVIRONMENT 

The antiferromagnetic system may be affected by a 
dissipative dynamics due to the interaction with its en- 
vironment. In principle the most common source of dis- 
sipation will be phonon excitations in the lattice. Thus, 



the interaction Hamiltonian will be given typically by the 
so-called spin-boson model [5, 18] 

V = EX>(<"i)( 5 ? + SJH + S* r ){A r , 3 + 4, .), (16) 



where A and are the annihilation and creation opera- 
tors of the environmental boson modes (which are consid- 
ered to be local) , and we have assumed that the coupling 
function g(ujj) is isotropic and the same for every mem- 
ber of the lattice. In addition, the Hamiltonian of the 
environment is 



3 r 

which is written as 

H E = J2I2 OJ i A L Ak >i> 



(17) 



(18) 



after taking Fourier transform. 

In linear spin wave theory approximation, the interac- 
tion term reads 



Vlsw= E J E^(^)[(l-i)VfK + ^) 

+ ( 1 + i )\/f (<4 + &r)+ 25 



-(<4a r - b\,b r )] (A rJ + Alj). 



Now the whole Hamiltonian has become much more in- 
volved than the original spin- wave theory Hamiltonian. 
However, we can consider a simplified version of this 
Hamiltonian based on the following facts, 

• We ignore the term 2S(A r j +A\,^) because it com- 
mutes with all observables of the system, which are 
the only interesting for us, as we going to trace out 
the environment later. 

• The terms a\,a r and b\.b r are negligible in compari- 
son to the others in the regime where the spin wave 
theory is valid (6). 

• Since |(1 ±i)/>/2| = 1, we have taken all the coef- 
ficients of a r , aj., 6 r , &J. to be just y/~S. This is not 
really a substantial modification to the dynamics 
but simplifies subsequent computations. 

Therefore, we arrive at 

Vlsw = y/s^^giwj) K + b t + 4 + K) (Arj+Alj), 

3 r 

(19) 

and by taking Fourier transform this interaction becomes 
V L sw = VS^j E fc 9^j) [(««, + &L)(^-fcj + 4,,) 



+ (a{ + b k )(A ktj + Al kJ ) 



(20) 
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Finally, after the Bogoliubov transformation of Eqs. (11) 
and (12) we get 

x ■(a fc + /3i)(A_ fc|J - + 4 |J .) (21) 



IV. MASTER EQUATION FOR A THERMAL 
ENVIRONMENT 



The dynamics of system and environment is given by 
the von Neumann equation 



where 



(22) 



(23) 



We aim at writing a dynamical equation for the state of 
the system ps = Ti E (p), where the trace is taken over the 
environment degrees of freedom. This task is generally 
quite complicated. However, we are particularly inter- 
ested in describing how the system evolves to the Gibbs 
state because the lack of insulation; and such a case is 



expected to happen for a large environment in thermal 
equilibrium (a "bath") with a small coupling constant. 
Under these conditions an equation, called master equa- 
tion, can be found by resorting to perturbation theory 
[19]. 

The initial state of the environment is then written as 



Pe 



Z-V^ Efe ^ A ^ Afe ^', (24) 



where Z — Tr (e @ He ) is the partition function with f3 = 
1/k^T. From now on we shall use natural units h = &b = 
1. 

Because of the Riemann-Lebesgue lemma [4] , for small 
coupling g(ujj) we can safely neglect the counter-rotating 
terms in (21), 



1/4 r 



V LS w= V^E.E,^)^) 



(25) 



Now the problem becomes equivalent to two collections 
of uncoupled harmonic oscillators given by their opera- 
tors a k and j3 kl which are coupled to set of independent 
environments characterized by k. The standard tools to 
obtain a master equation for a weak interaction with the 
environment can be found in references [4-7] . If we apply 
those techniques to this system we arrive at 



~j t = £(p) = ~ i[H LSW ,p] 

+ ^2lk(n k + 1) (a k pa{ - ]:{a ] k a k , p} + (3 k pp{ - Ufa*, p}) 
k ^ ' 

+ lkfi k {a{pa k - ^{a k a{, p} + p{p/3 k - ^ pj^j 

+ lk{n k + 1) (a k pp ] _ k - ^{/3l k a k , p} + p- k pa ] k - i{a£/?_ fc , p}^j 

+ lkfik (fi- k pOL k - ^{a fc /3i fc , p} + a ] k pl3- k - ^{f3- k a{, p}^j . (26) 

I 



Here 



7fc:=2 ^VdTf JMfc)) ' (27) 

where J(uo) = ^ ■ g 2 (u)S(uj — ujj) is the so-called "spec- 
tral density" of the bath. This one, for solid state en- 
vironments, is usually parameterized in the continuous 
limit [5, 18] as 

J(u) = auj s uj s c - 1 e- Uj/uJc , (28) 



where a accounts for the strength of the coupling and uj c 
is the cut-off frequency of the bath. Typically three cases 
are distinguished, s > 1 (super-ohmic), 8 = 1 (ohmic) 
and s < 1 (sub-ohmic) . The other quantity n is the mean 
number of phonons in the bath with frequency u;(fc), 

n k := [eMu(k)/T)-l}-\ (29) 

Note the presence of the crossing terms "a(-)/?t" in the 
master equation. They are absent in the analogous treat- 
ment of a ferromagnetic system and they will be induce 
a different rate for the decay of the magnetic order, for 
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k z = k z = ±7r 




FIG. 2: Magnon decay rate (27) in the first Brillouin zone. 
On the left it is depicted the surface for k z = and k z = ±7r 
is on the right. 



instance. 



A. Approach to the Equilibrium 

By construction [19], the Gibbs state p t h = 
Z~ 1 e~ HLSW / T , at the same temperature T as the bath, 
is the steady state of equation (26), i.e. £(pth) = 0. This 
is straightforwardly verified by taking into account that 

e -H L sw/T ak = e "W/T ake -H LSW /T^ (3Q) 
e -HL8w/Tp k = e «,(h)/Tp ke -H L sw/T m (31) 

Moreover any initial state of the system becomes closer 
and closer to this Gibbs state during time evolution. 

We have thus constructed a dynamical equation to de- 
scribe the thermal relaxation process of a quantum an- 
tiferromagnet. Remember that for the spin wave theory 
to make sense the number of magnons has to be small 
(6), so for large bath temperatures this treatment is not 
valid in the long time limit where the state approaches to 
the Gibbs' (which contains a large number of magnons 
for large T). However the predictions of equation (26) 
should also agree reasonably well with the exact ones at 
short times. 



B. Magnon Decay Rates 

A remarkable property of this system is that every ex- 
ponent is not allowed in the spectral density (28) in order 
to obtain finite results for many-body observables. This 
is because quantities such as magnon decay rates (27) 
and the thermal number of phonons become infinite for 
certain values of fc, so for those values the spectral den- 
sity has to approach zero fast enough. Particularly, it 
requires a super-Ohmic spectral density. It is worth to 
remind here that these kind of problems may also arise 



in simpler systems, for instance in a single spin when 
subject to a pure dephasing environment (see [4]). The 
concrete values of the rest of parameters of J(uu) is not 
very relevant for our proposes as we always assume to 
be in a sufficiently weak interaction regime [20], we shall 
take 

s = 3. a = J/10, u c — maxcj(fc) = 2JSd. (32) 

k 

In figure 2 we have represented two sheets of the 
magnon decay rate (27) in the first Brillouin zone. On 
one hand, we note that the magnon decay rate vanishes 
on the origin and on the eight corners of the Brillouin 
zone k — (±7r, ±7r, ±7r). On the other hand it reaches 
the maximum value on the points of a sphere of ra- 
dious r ~ 0.947 centered just this eight minimum points 
k = (d=7r, ±7r, ±7r). Between both cases there is a transi- 
tion which we have tried to illustrate by taken the values 
k z = 0, ±7r in the figure (given the symmetry of the decay 
rate, we can use k XjV instead of k z leading to the same 
figures). 

From (26) it is possible to compute the evolution of 
any combination of and in the Heiserberg picture. 
We give the result for the evolution of the first and sec- 
ond moments in Appendix A. Those expressions allows 
us to compute the evolution of any spin operator in the 
quantum antiferromagnet and relevant observables con- 
structed out of them. 



V. DYNAMICS OF RELEVANT OBSERVABLES 

In this section we study the time-evolution of some 
properties which have special interest in the description 
a quantum antiferromagnet. For concreteness, we have 
selected two of them, the staggered magnetization and 
the spin correlation functions. 

A. Staggered magnetization 

Due to the isotropy of Hamiltonian (1), one may expect 
the ground state to be also symmetric under rotations. 
However, as we have already mentioned, the ground state 
turns out to be close to the Neel state, which has clearly 
a privileged orientation. This is an example of sponta- 
neous symmetry breaking [14] . The figure of merit to 
compute this order in a quantum antiferromagnet is the 
expectation value of the staggered magnetization opera- 
tor 

™? = ^£(-l) IHI ^, (33) 

r 

which in the thermodynamic limit reads 

m st = lim (ml) = lim 1 £(-l)IMI<S*>. (34) 
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By using the equations (2) and (3), we may write this 
operator as 



with hk 



N 

r 

s - ^- n r = s — -J- 

N ^ N ^ 

r k 



(35) 



n^ ; . Note that k = 27rm/(N/2) where 



m varies two by two instead of one by one. So in the 
thermodynamic limit 
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Time (1/J) 



^Jim — y ^ = 2/2^)3 / (^ft) FIG. 3: Decay of the staggered magnetization showing the 

°° k ^ / ./B.Z. approach to the Gibbs state values mf. The red line cor- 



Here B.Z. denotes the first Brillouin zone, and the extra 
factor 1/2 appears because the double spacing between 
consecutive k on the left hand side. Thus, the staggered 
magnetization becomes 



16tt 3 



dk(h k ). 



(37) 



B.Z. 



When the system is interacting with a thermal bath, 
the staggered magnetization approaches in time to its 
thermal value. This is exactly zero for any T =^ due to 
de Mermin- Wagner theorem [21] in ID and 2D, however 
that is not the case in 3D. Additionally, note that the 
interaction Hamiltonian (16) is also isotropic so it is not 
trivial to find also magnetic order when the quantum 
antiferromagnetic is not isolated. 

From the master equation (26) we are able to visualize 
how staggered magnetization varies as a function of time. 
For this aim we just need to find the evolution of the 
observables hk- In terms of the operators ak and we 
have 



cosh 2 (6 k )ala k 

sinh(^)cosh(^)(4^ 

a k pk) + sinh 2 (# fc )(/^ fc + 1), (38) 

s'mh 2 (6 k)(ala k + 1) 

sinh(#fc) cosh(0 fc )(4/j£ + ^kPk) 

cosh 2 (#fc)/^ fc . (39) 



Particularly, if we start from the ground state, 
(a{a k (0)) = (4/3 fe (0)) = (a k p h (0)) = 0, we find 



n (6) 



(n { k a \t)) = (ti£>(t)) = cosh(20 k )fik [l - e 
+ sinh 2 (<9fc). 

Introducing these values in (37) and using equation (13): 



-iht . 



3(7**)] 



St 



m; 



St 



d 

"8^ 



f dk 



-iht 



cos(7fc£)] , 
(40) 



responds to T - 
T — IK and m s , 
mf ~ 0.237. 
short times. 



0.9 K with mf ~ 0.302, for the green line 
~ 0.271, and for the blue line T = 1.1 K and 
The inset shows more in detail the evolution at 



where 



= S 



16tt 3 



L 



dk 



(41) 



is the expectation value of the staggered magnetization in 
the ground state. In 3D, for a square lattice and S =1/2 
this value is rajf ~ 0.422. 

In figure 3, the evolution of the staggered magnetiza- 
tion is shown for different values of the bath temperature. 
It is noteworthy to mention the non-exponential decay of 
m st . This is due to its dependence on t through the inte- 
gral of (40), which renders combinations of different ex- 
ponentials. Remarkably, there is a short period where the 
order is lost very fast (between t = and t ~ 0.1/J, see 
inset figure). After that, the system continues evolving 
slower to the Gibbs state. This suggests that if we want 
to visualize variations of m st due to the environment, the 
best chance is to look for them in systems with not very 
small J. 



B. Two-point correlation functions 

It is also worthwhile to study the second moments of 
angular momentum operators. For the sake of illustra- 
tion, we focus in this section in the transversal two-point 
spatial correlation function, which is 



S ± (r 1 ,r 2 ,t) = Tr[SZ 1 SZ 2 p(t)}. 



(42) 



Without lost of generality we take r\ G A. Then for 
r 2 G A 

S 

S±(r 1 ,r 2 ,t) = -Tr[(a ri +4 1 )(o ra +a+ 2 )p(t)], (43) 



and 



S±(r 1 ,r 2 ,t) = -Ti[(a ri +4 1 )(6 ra + &+>(*)], (44) 
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for r 2 € B. 

Details of the computation are found in Appendix B, 



finally in the thermodynamic limit we obtain 



S±(ri,r 2 ,t) = 



S 



2(2tt) : 



/ dk cos [fe • 
Jb.z. 



(n - r 2 ) 



6(r 2 ){l + 2n fc [l - e-T fct cos( 7fc t)]} + 6>(r 2 )2n fc e"^* sin( 7fc t) 




= 
- - -t = 10" 3 

-♦-t = io- 2 

i = lO" 1 




-10 -5 5 



• T = K 

- - -T= 1 K 

- ♦ -T = 5 K 



-10 -8 -6 



-2 2 

n - r 2 



6 8 10 



FIG. 4: Evolution of S±(n,r 2 ,t) for T = 5 K in the ther- 
modynamic limit for different time instants (in units of J -1 ). 
Note the oscillating behaviour typical of antiferromagnetics 
systems. The inset figure illustrates the similarity between 
the cases for \S±(r±, r2, t)\ after normalization. 



FIG. 5: Thermal values of S±(ri, 7*2). In the inset is repre- 
sented again its normalized absolute value. 



the z direction the inelastic scattering is related to the 
correlation (S^_ k (t + T)Sjj,(t)) , where 



where 



and 



e fc (r) = 



e fc (r) 



d, if r e A, 
if r e B, 



if r e A, 
d, if r G 5. 



(46) 



(47) 



We have plotted this correlation for some time instants 
in figure 4. In addition, figure 5 shows different cases 
when the Gibbs states has been reached. 



C. Response function 

Other interesting quantities in this system are the re- 
sponse functions. They are the Fourier transform of two- 
time correlation functions of spin operators, and for in- 
stance, they directly appear in cross sections of inelastic 
neutron scattering, which are experimentally accessible. 
For an antiferromagnet with staggered magnetization in 



QX 



1 



ikr 



QX 



(48) 



One has to be specially careful when computing 
multitime-correlation functions for non-unitary evolu- 
tions. This is because the evolution of the product of two 
operators, say a and 6, does not equal to the product of 
the individual evolutions of a and b when the dynamics 
is not unitary, i.e. (ab)(t) 7^ a(t)b(t). However we can 
circumvent this problem by writing the correlation func- 
tion on the extended space where the evolution is indeed 
unitary 



\si h (t+T)sm\o) 



Tr e iff( * +T) S* fc e 



Tr([S*_ k (t 

Ht 



T)sz(t)\o)(o\® PE ] 

Ht\ 



] PE 



Here the trace operation is taken over both system and 
environment degrees of freedom, and H = Hlsw +He + 
Vlsw is the whole Hamiltonian of system and environ- 
ment. Then it is possible to obtain that (see detailed 
discussion in [7]) 

(0\Sl k (t + r)Sm\0) = <0| [S x - k (r)S x k ] (t)|0>. (49) 
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That is, it is needed to obtain first the Heisenberg evolu- 
tion with respect to the parameter r of the operator S__ k 
and after that, the Heisenberg evolution with respect to 
the parameter t of the product S__ k {r)S k . 
For linear spin-wave theory we have 




if r e A, 
if r e B, 



(50) 



thus, according to (7) and (48), 

S x k = ^(a k + al k + b_ k + b{). (51) 

If we perform the Bogoliubov transformation (11) and 
(12), the Eqs. (Al) and (A2) lead to 



S x - k (r) = ^e-^/ 2 [cosh(0 fc ) - sinh(0 fc )] (V^*^ {[cos( 7fe r/2) + sin( 7fc r/2)]a_ fe + [cos( 7fc r/2) - sin( 7fc T/2)]/3 fc } 
+ e [ ^ T {[cos( 7fc r/2) + sin( 7fe r/2)]4 + [cos( 7fc r/2) - sin( 7fe r/2)]/3l fc }) , (52) 

where we have used the fact that 6L& = k and uj(—k) = uo(k). Since by assumption j k is small, for small r we can 
neglect it in comparison to the complex exponential e ±1UJ ^ r , 



S x - k (r) 



cosh(^fc) — sinh(^fc)] 



-(a_ k +(3 k ) +e i ^(4 + /3l fc ) 



(53) 



Finally, by using (A4) and (A5), we compute the evolution of S__ k (r)S k with respect to t, and after simplifying 
vanishing terms the correlation function reads 

(S x - k (t + T)S x k (t)) = f [cosh(20 fc ) - sinh(20 fe )] {e^ ^[(a. k al k (t)) + + (al k (3 k (t)) + <a_ fc /?J(t)>] 

+ e^[(ala k (t)) + (tf_ k P- k {t)) + (a{/3. k (t)) + (a fe /?l fe (*)>]} 
- S(d ~ a) (e--( fe )^ {n k [1 - e-T** cos( 7fc t) + 2 sm( lk t)] + 1 } 



e i "( fc )-n fe [1 - e-T fct cos( 7fc i) + 2sin( 7fc t)]) . 



(54) 



The Fourier transform with respect to r leads to the re- 
sponse function: 

S ± (k,t,w) = SZ SQ (k,t)5[u-w(k)]+S+ SQ (k,t)5[w+u>(k)], 

(55) 



with 
S 



LSQ 



(k,t) 



g(d-gfc) 



+ 



-2 {n k [1 

2sin( 7fc i)] + l}, 
S(d-£ h ) 



cos(7 fc t) 



(56) 



n k [l - e 7fc£ cos(j k t) 



+ 2siii(7 fc t)]. (57) 

Therefore, at t = (or T = 0) only the form factor 
S^ S Q(k^t) remains. On the other hand we conclude 

that as temperature increases, S^ S g(k^t) also increases, 
and they have the same geometry in momentum space 
as the magnon decay rates 7&. However the behaviour 
with time is not monotonic, S starts increas- 

ing with time until it reaches its maximum value at 
t max = arccos(l/V^10)/7fc, after this point it decreases 
monotonically towards the thermal value. 



VI. CONCLUSIONS 

In this work we have analyzed the behaviour of a quan- 
tum antiferromagnet in contact with a boson thermal 
bath. Based on spin wave theory, we have applied the 
weak coupling procedure (Davies' theory) to obtain a 
master equation for the dynamics. We belief this is a ba- 
sic and fundamental problem which has remained quite 
unexplored so far. It is at the crossroads of strongly cor- 
related systems and the physics of quantum open systems 
that is so much rooted in quantum information theory. 

From the open systems point of view, spin wave the- 
ory provides us with a nice framework to apply the well- 
known techniques developed for quantum optics or quan- 
tum chemistry settings to quantum many-body prob- 
lems. Interestingly, some features which are typically 
encountered in small systems under weak coupling limit, 
e.g. the exponential decay of observables, may be lost 
when computing the observables which are relevant for 
the many-body systems. We have exemplified this point 
by studying the staggered magnetization, which for mod- 
erate temperatures, and despite of the the isotropic cou- 
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pling to the bath does not vanish. In fact, it does not 
shows an exponential decay either. 

Furthermore, we have illustrated the versatility of our 
master equation approach to the dynamics of thermal 
effects in quantum antiferromagnets by computing two- 
point correlation and response functions, also known as 
form factors. The geometry in momentum space of these 
response functions SLSQ(k,t) are closely related to that 
of the decay rate function in the first Brillouin zone. 
These form factors, in turn, are directly related to differ- 



current knowledge of a quantum antiferromagnet under 
non-isolated situations. 
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Appendix A: Time-evolution of the first and second moments 

The generator in Heisenberg picture is obtained by the equality Tr[X C(p)] = Tr[p£$(X)], for any operator X. 
By solving the dynamical equations, we obtain 



a k {t) 
<*kPl k (t) 



-ia;(fe)-7 fc /2]t 



cos(7 fc t/2)a fc (0) - sin( 7fc */2)/?_ fc (0)}, 



e [-i.( fc )- 7fc /2]t [sin(7 ^ /2)a _ fc(0) + cos ( 7fct/2 )^(0)], 
\e~^ {[cos( 7 ^) + l]a fc /?I fc (0) + [cos( 7fc *) - l]4/?_ fc (0) 
in( 7fc *)[a£a fc (0) + ^l fc /3_ fc (0) - 2n k }} , 
\e-^ 1 {[cos( 7fc t) + l]a{a k (0) + [cos( 7fc *) - l]^ fc /3_ fc (0) 
hkt)[a k (3l fc (0) + 4i8-fc(0)]} + n k [l - e"^* cos( 7fc *)], 
{[cos( 7fc t) - l]al k a_ k (0) + [cos( 7fc t) + l]^/3 fc (0) 

+ sin( 7fc *)[a_ fc /?£(0) + al fc /? fc (0)]} + n fc [l - e"^* cos( 7fc t)], 

a fc a fc (f) = | e-^W+^J* [(e^ £ -l) 2 /3_ fc /3_ fc (0) + (e^ +1) 2 <W0) - 2(e 2 ^ £ - l)a„/3_ fc (0)] , 

/W*) = | e-^W+^J* [(e^ -l) 2 a_ fc a_ fc (0) + (e^ +1) 2 /W0) - 2(e 2 ^ - l)a_ fc /3 fc (0)] , 

a k /3- k (t) = \ e -2(iwW+7i.)t [ (e 2 7fc t _i Kafc(0 ) + ( e 2 ^ -l)/3_ fc /3_ fc (0) - 2(e 2 ^ +l)a fc /3_ fc (0)] , 



tow 



2 K 

sin( 

§< 

sin l 



2 e 



-27fct 



(e^* -l) 2 /3_ fc /3+ (0) + (e^< +l) 2 a fe a t _ fe (0) - 2(e 2 ^' -l)a*4(0) 



i8_ fc ^(t) = | e" 2 ™* [(e™* -l) 2 a fc al fc (0) + (e^< +l) 2 /3_ fe /5+ (0) - 2(e 2 ^< -l)a fc 0j(O) 



<7 fc a_ fc (i) 



(e"27 fc t _i) afea t_ fc(0) + (e -2 7fct _i)^_ fc ^t (0) + (e -2 7fe * + i )afc/3 t (0) ' 



(Al) 
(A2) 

(A3) 
(A4) 

(A5) 

(A6) 
(A7) 
(A8) 
(A9) 

(A10) 

(All) 



1 -2(iw(fe)+7„)t \( P 1kt 



[(eT** +l) 2 a fc a_ fc (0) + (e™* -l) 2 /3 fc /?_ fc (0) - (e 2 ™* -l)(a fc /3 fc + a_ fc /3_ fc )(0)] (Al2) 



1 -2(iw(fc)+ 7fc )t [/ 7fc t 



[(eT** +l) 2 /? fc /?_ fc (0) + (eT** -l) 2 a fc a_ fc (0) - (e 2 ™* -l)(a fc /3 h + a_ fc /3_ fc )(0)] (Al3) 



= | e- 2 ^^)' [(e^ +l) 2 a fc /? fc (0) + (e^ -l) 2 a_ fc /?_ fc (0) - (e 2 ^ -l)(a fc a_ fc + /5 fe /3_ fe )(0)] (A14) 



Appendix B: Computation of the 2-spin correlator S±(ri,r2,£) 
We can compute the evolution of S±(r 1^2) in the Heisenberg picture. For a ri a]. we have 



(a ri ai 2 )(t) = ^£e- ifcr V fc '^aU*) 



(Bl) 



k,k' 



L J2e- ik ^e ik '^{[cosh(0 k )a k - sinh(^)^][cosh(^)4' - sinh(0 fc ,)/MX*)- 



k,k' 
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By using Eqs. (A1)-(A5) and simplifying the mean values which vanish on the ground state, the sum is left only with 
the terms where k = k': 

(a ri 4 2 (t)) = ^Y, e ~ ik<ri ~ r2) l cosh2 (°^( a i a ^ +s ^ 

k 

= l^Y / e ~ ik ' iri ~ r2) {cosh(20 fc )n fc [1 - e-T**cos(7fct)] + cosh 2 (0 fc )} . (B2) 

k 

Similarly, the remaining terms for r 2 G A give: 

= j^J2 e ~ ik ' ir2 ~ ri) l cosh2 ^( a k a *(t)) + sinh 2 (^)(/3t/3 fc (t)) + S inh 2 (^)] 

A k 

= _L ^e-*^" 1 - 1 ) {cosh(2^)n fc [l - e"^ cos( 7fc t)] + sinh 2 (^)} , (B3) 

k 

(a ri a r2 (t)) = (at i at 2 (t))* = ^^e- ifc -^-^)sinh(2^)n fe e-^ S iii(7 fe i). (B4) 
and for r2 G 5: 

(a ri b r2 (t)) = (atXM* = E e " ifc ' (r2 " ri) cosh(^)sinh(^)[(4a fc (t)) + (fifo®) + 1] 

= -7=i=^e- ifc -^-^)sinh(2^){2n fc [l - e"^* cos( 7fc t)] + 1} , (B5) 
<a ri &t 2 (*)> = = ^== J2e~ ik < r ^ cosh(26» fc )n fc e-^sin( 7fc t). (B6) 

Thus, on one hand, the correlation function for r 2 ^ ^ reads 

S±(ri,r 2 ,t) = ^cos[/c • (n - r 2 )] (cosh(20 fc ) {2n fc [l - cos{ lk t)] + 1} - 2 sinh(20 fc )n fc e^ fct sin( 7fc *)) 

+ i^-E sin [ fc -( r i-^)]- (B7) 



Since 



^- ^sin[fe • (n - r 2 )] = Im |^ E eifc ' (ri " r2) } = M<W = °> ( B8 ) 

and using (13): 

5_L(ri, r 2 , t) = 2^ cos[fc • (n - r 2 )J I ^ rf2 _ ^ 2 I • (B9) 

And on the other hand, for r 2 € -B, 

q / -5 r . , + 2n fc {&[! - e~^ cos( 7fc t)] - rie~^ sin( 7fc rQ} ^ 

5±(ri ' r2 '* ) = ^wrm^ cos[k ' {ri ~ r2)] { T^i ) ■ (B10) 
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